Study of Unsteady, Sphere-Driven, Shock-Induced Combustion 
for Application to Hypervelocity Airbreathing Propulsion 
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Abstract 

A premixed, shock-induced combustion engine has been proposed in the past as a viable 
option for operating in the Mach 10 to 15 range in a single stage to orbit vehicle. In this 
approach, a shock is used to initiate combustion in a premixed fuel/air mixture. Apparent 
advantages over a conventional scramjet engine include a shorter combustor that, in turn, 
results in reduced weight and heating loads. There are a number of technical challenges that 
must be understood and resolved for a practical system: premixing of fuel and air upstream 
of the combustor without premature combustion, understanding and control of instabilities 
of the shock-induced combustion front, ability to produce sufficient thrust, and the ability to 
operate over a range of Mach numbers. This study evaluated the stability of the shock- 
induced combustion front in a model problem of a sphere traveling in a fuel/air mixture at 
high Mach numbers. A new, rapid analysis method was developed and applied to study such 
flows. In this method the axisymmetric, body-centric Navier-Stokes equations were 
expanded about the stagnation streamline of a sphere using the local similarity hypothesis in 
order to reduce the axisymmetric equations to a quasi-ID set of equations. These reduced 
sets of equations were solved in the stagnation region for a number of flow conditions in a 
premixed, hydrogen/air mixture. Predictions from the quasi-ID analysis showed very 
similar stable or unstable behavior of the shock-induced combustion front as compared to 
experimental studies and higher-fidelity computational results. This rapid analysis tool 
could be used in parametric studies to investigate effects of fuel rich/lean mixtures, non- 
uniformity in mixing, contaminants in the mixture, and different chemistry models. 

Nomenclature 

= Local speed of sound 
= Forward reaction rate coefficient 
= Temperature exponent for forward reaction rate 
= Damping term 
= Total energy 
= Activation energy 
= Total enthalpy 
= Sensible enthalpy 
= Flux vectors 
= Source vector 
= Sphere radius 

= Reynolds number based on freestream quantities 
= State vector 
= Mass fraction 
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k = Local body curvature 

M = Mach number 

n, s = Normal and tangential coordinates 
p = Pressure 

q = Heat flux 

r = Local cylindrical radius 

R u = Universal gas constant 

t = Time 

T = Temperature 

v, u = Normal and tangential velocity components 
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Body angle measured from stagnation streamline 
Non-dimensional height measured from stagnation streamline 
Ratio of specific heats 

Non-dimensional distance measured from sphere center 

Dynamic viscosity 

Production rate of species i 

Stress tensor 

Density 


I. Introduction 

Access to higher Mach numbers is one of the goals that drive current research in supersonic combustion ramjets 
(i.e scramjets). The benefit of operating at higher Mach numbers is evident when considering the example of a 
scramjet that is integrated in a single stage to orbit space-access vehicle 1 . For example, at Mach 10 such a vehicle 
has attained approximately 15% of the required kinetic energy to reach low Earth orbit (LEO) before transition to 
rocket-mode while at Mach 15 it has reached 30%. To date, most research and design concepts are projected to 
operate up to Mach 10 before experiencing a significant reduction in or loss of performance. This is because 
increasing the freestream Mach number reduces the time that is allowed for diffusive mixing and combustion. One 
of the simplest ways to mitigate this problem is by elongating the combustor; however, this increases the mass of the 
vehicle considerably as the combustor tends to be a significant mass element of the engine. 

The utilization of premixed, shock-induced combustion (PMSIC) is one method that shows promise in driving 
airbreathing propulsion to higher Mach numbers (i.e. hypervelocity propulsion). In a PMSIC engine concept (see 
Figure 1), the fuel is injected on the forebody of the vehicle rather than in the combustor, as is the case for 
conventional scramjets. The fuel and air undergo mixing and then travel through a strong, oblique shock wave from 
the cowl lip that causes auto-ignition and subsequent combustion of the fuel-air mixture. Apparent advantages of 
the PMSIC engine concept are a shorter combustor length that, in turn, results in reduced weight and heating loads. 
Assuming the vehicle is boosted to hypervelocity speeds, the isolator component typically found in scramjet 
concepts is able to be removed, further amplifying these benefits. Other benefits include improved pressure 
recovery at the combustor when hydrogen fuel is used 2 . These benefits translate into access to much higher Mach 
numbers than conventional scramjets and, therefore, delayed staging when used for space access. 

When considering the problem of injecting fuel in the inlet, there are multiple technical challenges that must be 
overcome to arrive at a practical concept. These challenges include 

• Inlet auto-ignition - Because the fuel is injected on the forebody of the vehicle, it is possible for the 
tiic l/air mixture to experience pre-ignition in the hot boundary layer on the surface. Large amounts of 
drag and heating are then experienced on the forebody, resulting in significant performance losses in the 
engine. 

• Mixing losses - The premixed fuel-air stream may not be uniformly mixed or may be lost due to inlet 
spillage. These non-uniformities may have an important effect on the stability of the ignition and 
combustion. 

• Shock wave stability - Depending on the flow conditions, a weak or strong coupling between the 
reaction front and the shock wave inducing combustion may occur, resulting in instabilities of the 
reaction front and the shock wave. The presence of these instabilities would have consequences 
ranging from reduced performance to engine unstart. 
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The focus of this paper is to examine shock wave stability in a representative model problem of hydrogen/air 
PMSIC ahead of a blunt projectile traveling at near-hypersonic speeds. Previous numerical and experimental work 
has been conducted on this particular problem with one observation being that instabilities in the flow emanate from 
the stagnation region of the body and propagate downstream 3 . This study takes advantage of this observation and 
re-formulates the blunt-body problem in a quasi-ID manner. In this formulation, the flow was numerically simulated 
only along the stagnation streamline with pertinent gradient values normal to the centerline included as well. This 
approach drastically reduced computational time over multi-dimensional simulations, thus allowing rapid studies of 
the effects of different parameters (e.g. free stream pressure, velocity, nose radius, varying stoichiometry, non- 
uniformity in mixing, contaminants, chemistry models, etc.). This rapid computational capability allows the analyst 
to quickly parameterize the behavior of the system over a wide range of design variables. This is important for 
identifying safe regimes where the flame may be stable for conditions found at the combustor during hypervelocity 
flight. As a result, this formulation should be a valid starting point for higher-fidelity flowpath and trajectory 
designs. For this study, the behavior of the reaction and shock fronts was evaluated over a wide range of velocities 
relevant to past experiments. The effect on stability of changing inflow stoichiometry was also studied. 

II. Body-Driven, Premixed Shock-Induced Combustion 

To date there have been many experimental and numerical studies of the blunt-body PMSIC phenomenon. Early 
experimental work in the 1960's and 70's 4 " 7 provided observations of PMSIC in a controlled environment and also 
classified different regimes of combustion and shock instability. The first regime is a stable regime that exhibits 
smooth, uncoupled reaction and shock fronts with no unsteady motion. The stable regime will not be discussed 
further in this section. The next regime is a regular, unsteady regime in which a "corrugated" reaction front appears 
that is periodic and forms a small-amplitude oscillation in the bow shock. Finally, a large -amplitude, unsteady 
regime is observed that is non-periodic and exhibits significant coupling of the bow shock and reaction front such 
that the reaction zone pushes the bow shock far upstream before retreating back to its original position. Examples of 
these regimes may be found in Figures 3, 4, and 5 and will be discussed in greater detail later. 

The mechanism of the regular, unsteady regime was initially synthesized from experiments by McVey and 
Toong 5 and was later refined by the numerical studies of Matsuo and Fujiwara s to include the effect of reflected 
compression waves from the body. In this regime, pressure waves created by the inception of a new reaction front 
travel both toward the body and the bow shock. Collision of the pressure wave with the bow shock pushes the shock 
upstream and reflects a contact discontinuity. Upstream of the contact discontinuity is gas with higher temperature 
(and therefore shorter ignition delay time) which eventually results in the formation of a second reaction zone ahead 
of the original one and the production of another set of pressure waves. Once the contact discontinuity meets the 
original reaction front it is extinguished, thereby generating expansion waves. This behavior is repeated periodically 
in a staggered manner. 

A large-amplitude, unsteady mechanism was proposed by Alpert and Toong 7 along with experiments 
documenting its behavior. This model was later updated through numerical studies of Matsuo, et al. 9 It was 
proposed that the instability exhibited in this regime finds root in the significantly higher heat release caused by the 
harsher free stream conditions. This heat release at the reaction inception point causes a superdetonation to travel 
towards the bow shock and a retonation to travel toward the body. In this context, a retonation is a strong pressure 
wave that travels downstream through the burned gases while a superdetonation occurs where the reaction front and 
its initiating shock wave are tightly coupled. Interaction of the superdetonation with the bow shock causes a 
retonation and contact discontinuity to travel toward the wall while the superdetonation heaves the bow shock 
upstream in a large-amplitude manner. Upon relaxation of the bow shock and the reaction front, a mechanism 
similar to that of the regular, unsteady regime is observed in the creation of a new reaction front. This cycle repeats 
non-periodically with an effective period greater than that of the regular regime. 

Other numerical studies have shed light on the physical behavior and computation of body-driven PMSIC. 
Wilson and MacCormack 10 found in simulating the experiments of Lehr' 1 that there was significant sensitivity both 
to grid spacing and to reaction rates in the Jachimowski model 11 . Another study by Wilson and Sussman 12 applied a 
logarithmic transformation to the species continuity equations to provide enhanced solution accuracy for coarser 
grids. Investigations by Ahuja, et al 13-15 simulated the Lehr experiments at and above the detonation velocity of the 
mixture (Mach 5.11 and 6.46, respectively). In their studies, the Matsuo regular regime model for the Mach 5.11 
case was verified and the frequencies in the simulated flow field, determined using a Fourier transform, compared 
favorably with experiment. For the overdrive case of Mach 6.46 it was shown that the period of instability was 
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reduced to such an extent that the flow field appears macroscopically stable and the reaction front had merged with 
the bow shock. These numerical results were also verified against experimental data. 

III. Simulation Methodology 


Governing Equations 

The axisymmetric Navier-Stokes equations and species conservation equations in the body-oriented coordinate 
system shown in Figure 2 can be written as 


dU dM dN 

~rr + — h — F Q - 0 

at os on 
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The state vector (U), flux vectors (M, N), and source term vector (Q) are expressed in non-dimensional form as 
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where A = 1 + nk, (3 = r + n cos a. In these equations, variables have been non-dimensionalized with respect to 
free stream velocity, temperature, density, and body radius. 

In order to arrive at the quasi- ID form of equations, a local similarity postulation 16 17 was used in the stagnation 
region. Under this approach, flow variables were expanded in a series about the axis of symmetry by 


p(s,n,t ) = p 1 (n, t) 

p(s,n,t) — p-i (ji, t) + p 2 (n, t) sin 2 a 

H(s,n,t ) = 

u (s, n, t) = % (n, t) sin a 

v(s,n, t) = v 1 (n, t) cos a 

(First Truncation) 


+ p 2 (n, t) sin 2 a + ••• 

+ p 3 (n, t) sin 4 a + ••• 

+ H 2 (n, t) sin 2 a + ••• 

+ u 2 (n, t) sin 3 a + ••• 
+v 2 (n, t) cos a sin 2 a + 


(6) 


Note that the retention of the second order term in pressure is warranted by the order of magnitude analysis. The 
expansions given in Equation (6) have had a small-angle approximation applied due to the result that s — a for a 
non-dimensionalized sphere. The solution variables in the first truncation have a simplified notation where terms in 
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the expansion equal to zero have been removed and the remaining solution variables and their derivatives have 
simply been given the names of the expanded state variable with subscripts. Here it-, is actually the first derivative 
of velocity with respect to s and p 2 is the second derivative of pressure with respect to s. All other subscripted 
variables in the first truncation are equal to their respective state variables. 

Retaining terms only to the first truncation shown in Equation (6) and substituting into Equations (3) - (5), the 
following set of equations is obtained for a sphere by equating the coefficients of like powers of a. Note that on the 
stagnation streamline a — 0. The governing equations are reduced to the form 


dV dN' 
dt + dn + 


Q' 


= o 


(7) 


where now the solution is a function only of time and distance normal to the surface. Note that the primes have been 
added to indicate the quasi- ID form of the vectors. The state vector has the form 
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The flux vector is split into contributing terms by 
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where each term represents the inviscid, viscous, and chemical contributions to the flux, respectively. These terms 
are given by 
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In a similar fashion the source term is split up as 
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An additional equation is required for the solution of p 2 , given by 


dPz = dpi | + 

dn dn 1 + n 


(16) 


which comes from the ^-momentum equation. This equation is decoupled from the Navier-Stokes equations and can 
be solved separately in a non-lagging manner. 


Numerical approach 

The governing equations are integrated over the convective time step using two methods. First, the convective 
and diffusive terms in the governing equations are solved using MacCormack’s predictor-corrector technique 18 . 
This scheme is second-order accurate in both time and space. As the current implementation employs shock 
capturing, numerical damping is required to control the resulting dispersive error occurring around areas of large 
gradient in the solution such as the shock front. The terms describing numerical damping at each grid point have the 
TVD-like form 19 


where 
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and where S<p and are damping constants that have values the analyst may choose — for this study they were both 
set to 0.3 after some manual optimization. The placeholder variable <p represents any variable used in the damping 
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expression, for example pressure, temperature, or species concentration. If multiple variables are used in the 
damping scheme simply add each to the PDE, which at spatial step i takes the form 


/diry p?(N '-0 01 -Z) 02 

l dt I . dn 

x / j 




+ (Q ')) = o 




( 21 ) 


The damping term should be updated for each step of the MacCormack method. Note that when differencing the 
state vector in Equation (17), the opposite type of differencing should be used than what is used in the current 
MacCormack step. For example, if forward differencing is currently being used in the MacCormack step, 
backwards differencing should be used for the damping term and vice versa. 

Once the convective terms have been integrated, changes in species composition at each grid point due to 
chemistry are then computed. In the current implementation, the chemical reactions are integrated over the 
convective time step using a quasi-steady-state solution method described by Mott, et al. 20 The chemical model 
used in this study is a reduced form of the hydrogen-air model developed by Jachimowski 1 '. It is an 8-species, 13- 
reaction model which includes the chemical species H, O, OH, H 2 0, H0 2 , H 2 , 0 2 , and N 2 . The species N 2 only acts 
as a diluent species unless otherwise stated. The species H 2 0 2 is excluded from the present mechanism due to its 
slow rate of formation in high -temperature combustion applications 21 . Table 1 shows the mechanism used in this 
study and the parameters used in the Arrhenius equation for each reaction. 


Table 1. An 8-species, 13-reaction chemical model for H 2 -air combustion. 11 



Reaction 

A 

b 

Fyiooo 

1 

0 2 + H 2 ^ OH + OH 

1.70E+13 

0.0 

48.00 

2 

OH + H2 ±5 H + H 2 0 

2.20E+13 

0.0 

5.15 

3 

H + 0 2 £5 OH + 0 

2.60E+13 

0.0 

16.8 

4 

0 + H 2 £5 H + OH 

1.80E+10 

1.0 

8.90 

5 

OH + OH £5 H 2 0 + 0 

6.30E+12 

0.0 

1.09 

6 

OH + H + M £5 H 2 0 + M 

2.20E+22 

-2.0 

0.00 

7 

0 + H + M±5 0H + M 

6.00E+16 

-0.6 

0.00 

8 

h + h + msh 2 + m 

6.40E+17 

-1.0 

0.00 

9 

H + 0 2 + M ±5 H0 2 + M 

2.10E+15 

0.0 

-1.00 

10 

H0 2 + H ±5 H 2 + 0 2 

1.30E+13 

0.0 

0.00 

11 

H0 2 + H ±5 OH + OH 

1.40E+14 

0.0 

1.08 

12 

H0 2 + 0 ±5 OH + 0 2 

1.50E+13 

0.0 

0.95 

13 

H0 2 + OH £5 H 2 0 + 0 2 

8.00E+12 

0.0 

0.00 


Species: H, O, OH, 0 2 , H 2 , H,0, H0 2 , N 2 


The listed forward rate coefficients satisfy the Arrhenius fonn kf = AT b e R u T and are in cgs units. 
Third body coefficients for each reaction are as follows: 

(6) H 2 0 = 6; (7) H,0 = 5; (8) H 2 0 = 6, H 2 = 2; (9) H 2 0 = 16, H 2 = 2 


The inflow boundary condition simply takes on freestream values of the variables due to the method being shock 
capturing. u u p x , and 7\ were set to one, v 1 was set to negative one, p 2 was set to zero, and p 1 was obtained from 
the non-dimensional equation of state. The initial species composition across the domain was that of the freestream. 
At the wall, v 1 was set to zero due to impermeability, p 1 was set using a zero gradient assumption, and % was set 
either to be zero (viscous, no slip) or equal the value contained in the nearest grid point (slip). Temperature was 
either set using an adiabatic assumption or set to be constant when considering cases of a particular wall 
temperature. Density at the wall was determined using the perfect gas relation at the wall. 

The initialization strategy was the same for each case being considered. The values of wall pressure and both 
velocities were set to change linearly from their freestream value to a guess of their wall value. Total enthalpy was 
assumed to be constant along the domain at the freestream value, from which initial conditions of temperature, 
sensible enthalpy, density were then determined. Species composition across the domain was also assumed to be 
constant at the inflow value across the solution domain. 

A 400-point, equally spaced grid was found to be sufficient in this study in order to resolve the fine details of the 
flow. A CFL number of 0.4 with global time stepping was chosen for the results presented in this paper. All 
simulations were conducted on a high performance workstation with code written in Fortran. Because of the 
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unsteady phenomena present in the flow, limit cycling in the residual calculation is typically present with the 
exception of cases with a stable reaction front. The simulations were considered to be converged when the transient 
features of the flow were eliminated. 


IV. Results and Comparison 

Qualitative Experimental and Numerical Comparison 

The experiments of Ruegg and Dorsey 4 were simulated using the quasi-ID code in order to qualitatively assess 
the suitability of the proposed analysis method. Ruegg and Dorsey conducted early experiments on PMSIC by 
firing spherical, 20 mm diameter projectiles into a stoichiometric hydrogen-air mixture in the vicinity of Mach 5 at 
pressures ranging from 0.1 to 0.5 atm. Over the range of experimental conditions investigated, three distinct 
regimes of behavior of the reaction/shock front system were observed: stable; regular, unsteady; and large- 
disturbance, unsteady. A qualitative comparison between the quasi- ID simulation results and experimental 
shadowgraphs at the conditions in Table 2 will now be made. The simulated results consist of an evolution of water 
mass fraction and pressure on the stagnation streamline as a function of non-dimensionalized time and distance (n) 
from the surface of the projectile. 

Table 2. Selected Ruegg and Dorsey experimental conditions for a 20 mm diameter, spherical projectile 
fired into p remixed, stoichiometric hydrogen/air 4 . 


Case 


Poo (atm) 

Too (K) 

Stability Regime 

1 

5 

0.10 

300 

Stable 

2 

4.9 

0.25 

300 

Regular, Unsteady 

3 

4.3 

0.50 

300 

Large Disturbance, Unsteady 


Figure 3 shows quasi-ID simulation and experimental results for the stable regime of oscillation. The simulated 
and experimental results indicate reaction and shock front locations that are invariant with time. This invariance is 
observed experimentally by smooth shock and reaction fronts. The reaction and shock fronts are separated by a 
relatively large induction distance with no direct interaction between the two beyond the shock wave raising the 
temperature of the hydrogen/air mixture above its autoignition limit. There is qualitative agreement between 
simulation and experiment for this case. 

The regular instability regime is reproduced in Figure 4. The experimental shadowgraph shows a corrugated 
pattern in the wake of the sphere caused by an oscillation in density due to the periodic creation of new reaction 
fronts in the stagnation region of the flow. Between the corrugated, burned region of gas and the shock wave lies a 
region of unburned gas. This is due to the presence of an observable induction distance between the bow shock and 
the reaction front. Lines of density variation in the unburned region are due to periodic contact discontinuities 
traveling from the bow shock to the beginning of each new reaction zone. The contact discontinuities are not visible 
in the simulated results in this figure because contour lines of pressure are shown (pressure being equal on both sides 
of a contact discontinuity). The contact discontinuities in the computations will be visible in later figures when the 
computed results are compared to established ID mechanisms. There is qualitative agreement between the 
experiment and the present computations because the physical mechanism observed in the simulated results match 
the flow features witnessed in the experimental shadowgraph. 

The large-disturbance instability regime is shown in Figure 5. Experimental results of the flow at the Case 3 
condition show a large-amplitude heaving of both the bow shock and the reaction front. Both the bow shock and 
reaction zone oscillate with a much lower frequency than that observed in the regular regime and apparently non- 
periodic features are visible in the shadowgraph. Although both the experiment and the present numerical results 
show large amplitude behavior, it is difficult to isolate other features to determine how well the quasi- ID simulation 
is performing. For further comparison, the present results are compared to a higher-fidelity, 2D simulation. Figure 6 
shows a zoomed-in view of a section of the present quasi- ID simulation compared to a segment of a simulation at 
the same conditions conducted by Matsuo and Fujii 22 . These images show contours of density as a function of time 
and distance along the stagnation streamline. The present results show a striking agreement with the high-fidelity 
numerical result and clearly show coinciding features of reaction shocks, retonations, contact discontinuities, and a 
superdetonation as the bow shock is penetrated by the reaction front. These features will be discussed in further 
detail later, but for now it appears that the quasi- ID results are capable of predicting the large-disturbance instability 
regime. 


8 


American Institute of Aeronautics and Astronautics 



Investigation of Reaction Front Oscillation Frequency 

The current analysis approach has been applied to simulate the ballistic experiments conducted by H.F. Lehr 1 in 
which 15 mm diameter projectiles were fired into a stationary, premixed, stoichiometric hydrogen/air mixture at 
near-hypersonic speeds. The ambient conditions of the fuel-air mixture had a static temperature and pressure of 292 
K and 0.42 atm, respectively. Before the results from the current study are discussed, some quick observations are 
made about Lehr’s experiments and the numerical simulations by Sussman 2 '. 

Lehr investigated the behavior of the reacting system over a wide range of velocities, 1685 to 2058 m/s. For 
each case, the frequency of oscillation of the reaction front was determined using shadowgraph images and 
considering the distance between successive striations in the wake of the sphere. Lehr found that as projectile 
velocity increased the frequency of oscillation increased monotonically. At high projectile velocities an overdriven 
detonation occurred which caused an apparently stable reaction front. 

The Chapman-Jouguet detonation velocity of a stoichiometric hydrogen-air mixture with initial conditions as 
that of the Lehr experiments was calculated to be 1957 m/s using the Chemical Equilibrium with Applications 
(CEA) program developed at NASA Glenn by McBride and Gordon" 4 2 \ This velocity is sensitive to species 
present — in this case only the species H, O, OH, H 2 0, H0 2 , H 2 , 0 2 , and N 2 were allowed as product species. The 
experimental value reported by Lehr was 2030 m/s, corresponding to the velocity past which oscillations in the 
reaction front were no longer visible. It should be noted, however, that other analyses have established that 
oscillations might occur above the Chapman-Jouguet detonation velocity. For example, the experiments of Alpert 
and Toong 7 found oscillations in the overdriven range of flight and modern numerical analyses have also 
substantiated this observation 13 " 15 23 . 

Sussman 23 simulated the axisymmetric flow around a 15 mm diameter sphere at the same ambient conditions and 
freestream velocities studied by Lehr using an upwind, TVD scheme. The oncoming flow is composed of a 
stoichiometric hydrogen-air mixture. In general, the numerical results of Sussman matched well with the 
experimental results with the exception of cases near the high velocity limit. A large discrepancy was observed at 
2058 m/s, which was postulated to be due to the method of interpreting the experimental frequency at that point. 
However, the experimental oscillation frequency at this point was successfully replicated in numerical work done 
around the same time by Ahuja 13 ’ 15 . At a projectile velocity of 2119 m/s, Lehr showed no oscillation although the 
numerical work by Ahuja showed oscillations occurring all the way to 2605 m/s; however, at this point they were of 
very small amplitude and located away from the stagnation streamline. Sussman was able to simulate oscillations 
for the 2119 m/s case but not for the 2605 m/s case. 

Figure 7 shows results from the quasi-ID simulations over the velocity range studied by Lehr and Sussman. Note 
that due to the reduced computing requirements of the quasi-ID simulation, a large number of cases are computed 
over the velocity range of the experiment versus previous numerical studies. In doing so, a better idea as to the 
general shape and limits of the parameter space is obtained. It is observed that, in general, the present results are in 
reasonable agreement with experimental results and higher-order simulations. Overall the quasi- ID simulations 
over-predict oscillation frequencies by 1 7%-34%. In all cases the overall behavior of the flow is predicted (e.g. 
regime of instability) well. 

Table 3. Oscillation frequencies from experiment and numerical investigations. 

Projectile Velocity (m/s) Lehr 6 (kHz) Sussman 2 (kHz) Present Results (kHz) 


1685 

148 

135 

200 

1804 

425 

420/74 

546/89 

1931 

712 

711/103 

849/126 

2029 

1040 

1038 

1282 

2058 

1960 

1180 

1351 

2119 

n/a 

1500 

1744 

2605 

n/a 

n/a 

n/a 


One of the benefits of having the ability to more densely populate the projectile velocity parameter space is that 
transition points from one regime of stability to the next become clearer and the general morphology of the space 
begins to focus. In addition, transition points may be more finely interrogated by the quasi-ID simulation in order 
to observe the changing behavior of the flow in order to arrive at mechanisms for transition. With a more densely 
populated parameter space, trends not observed in previous work also appear. For example, the space takes on 
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something akin to a terraced appearance where there are distinct regions of stability with rapid transitions in some 
locations from one dominant frequency to the next (e.g. around 2000 m/s and 2080 m/s). 

The basic transition mechanism at work in the observed results is an out-of-phase behavior of pressure waves 
generated by the new reaction fronts. For projectile velocities below 1650 m/s, the reaction front is observed to be 
smooth. If new reaction fronts are being created, they are imperceptibly low in amplitude and strength at the studied 
grid resolution. At the beginning of the regular regime zone, there are moderate amplitude reaction front features. 
Increasing the projectile velocity has the effect of creating additional reaction zones along the same reaction 
peninsula. This is due to the pressure waves generated by new reaction zones being significantly out-of-phase with 
the reflected pressure waves from the projectile surface. As velocity increases, new reaction fronts are created by 
out-of-phase waves propagating faster than expansion waves can relax the shock front. This causes a secondary, 
low-frequency oscillation to appear mainly between 1730 m/s and 1960 m/s upon which the higher-frequency 
regular instability is superimposed. The secondary oscillation dies out around the detonation velocity of the mixture 
(1957 m/s), but this may be coincidental. As the projectile velocity increases past 1960 m/s, where the classic, in- 
phase regular regime is re-established, a “terraced” trend of oscillation frequency versus projectile velocity occurs. 
The cause of this trend is again due to narrow velocity ranges where out-of-phase pressure wave behavior creates 
new reaction fronts in addition to ones already being created periodically due to the classic regular regime 
mechanism. The reaction front at the transition points shows mixed regular regime behavior with non-cohesive, 
irregular reaction fronts created by the out-of-phase pressure waves. In general, as the projectile velocity increases, 
the amplitude of the reaction front is found to decrease, eventually smoothing out into apparently stable behavior as 
it begins to couple with the shock front. 

Comparison of Present Results to ID Mechanism Models 

The simulated experiments of Ruegg and Dorsey are now reexamined in order to compare the present quasi- ID 
numerical results to published ID instability mechanism models in the literature. First, the quasi-ID simulation of 
the regular, unsteady regime is compared with the published mechanism of Matsuo 2 ' 1 in Figure 8. Note that all of the 
major features of the mechanism proposed by Matsuo are present in the simulated result. A newly created reaction 
front produces pressure waves that travel downstream toward the projectile surface and upstream toward the bow 
shock. The upstream traveling pressure wave meets the bow shock roughly in phase with the pressure wave 
reflected from the projectile surface. The collision of the pressure waves and the shock front causes the shock to 
move upstream with a contact discontinuity reflected downstream. The contact discontinuity is visible on the 
density plot and is visible as a line traced from each cusp of the bow shock to the newly created reaction front. The 
hot gas upstream of the contact discontinuity has a lower induction time than the downstream gas, which creates a 
new reaction zone that eventually extinguishes the old one. This process repeats indefinitely for a projectile with 
constant velocity. The simulated result differs slightly from the published mechanism in that the simulated contact 
discontinuity is connected to the tip of the reaction zone peninsula versus the ID mechanism, which traces it to the 
point where the old reaction front is merged with the new one. In addition, the presently simulated result indicates 
that the reflected pressure waves skip a generation of new reaction zone before interacting with the bow shock. The 
published mechanism shows the reflected pressure wave interacting with the bow shock along with the upstream- 
traveling pressure wave of the next reaction zone generation. Finally, the present computations show a rarefaction 
wave reflected at the point where the pressure wave and the bow shock meet, which is not depicted in the published 
model. 

The quasi- ID simulation of the large-disturbance, unsteady regime, shown with a density contour plot, will now 
be compared with the corresponding mechanism of Matsuo and Fujii 2 " in Figure 9. The large-disturbance 
mechanism of Matsuo and Fujii requires an energy release large enough to create a deflagration to detonation 
transition (DDT), noted in the figure as an “explosion in an explosion” (this is discussed further in the detonation 
literature 27 ' 2S . This causes a penetration of the self-propelled reaction front through the bow shock, reflecting a 
rarefaction wave and contact discontinuity downstream at the point of intersection. The self-propelled reaction 
front, which is coupled with the bow shock for a short time, gradually separates from the shock wave and relaxes to 
its original position before resuming the cycle. Comparing the present results to the ID mechanism shows that all of 
the essential features of the mechanism are reproduced in the quasi-ID simulation. Figure 10 shows a detailed view 
of the detonation region of the quasi-ID simulation to examine what is happening closer to the shock wave 
penetration. Flere it is clearly visible that the intersection of the reaction shock wave from the first explosion with 
the bow shock is responsible for generating the rarefaction wave and the contact discontinuity that causes the second 
explosion. The second pressure wave seen in the simulated results is due to the generation of the second explosion, 
explaining why it doesn’t appear in the general, single-explosion mechanism. After the second explosion, the 
pressure wave and reaction front traveling upstream merge with the bow shock and a self-sustained superdetonation 
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is observed for a short amount of time. Overall the quasi-ID simulation shows excellent agreement with the ID 
mechanism found in the literature. 

Further Parametric Assessment 

While the quasi- ID formulation is capable of rapidly evaluating the effect of many parameters relevant to the 
stability of the reaction front, for brevity only one case study is shown here. A point of comparison with the 
previously presented Lehr results can be made by adjusting the freestream stoichiometry to be either lean or rich to 
study the effect on the behavior of the reaction front. For this comparison three cases were chosen across the 
velocity range studied by Lehr and the freestream stoichiometry was set to 0.80 for the lean case and 1.20 for the 
rich case. All other quantities (e.g. velocity, pressure, temperature) were held the same as in experiments. In 
general, for hydrogen-air, it was found that lean freestream conditions shifted the behavior of a particular projectile 
velocity case to a higher effective velocity while rich conditions shifted the behavior to a lower effective velocity. 
For example, for a projectile velocity of 1696 m/s, the stoichiometric freestream behavior is in a regular unsteady 
regime. Shifting the freestream stoichiometry to a value of 0.80 results in two frequency behavior while shifting the 
stoichiometry to 1.2 causes the reaction front to be completely stable. This behavior is due to the fact that the 
molecular weight of a lean fuel-air mixture increases for a lighter-than-air fuel such as hydrogen and the gas 
constant of the mixture increases accordingly. Because temperature, pressure, and velocity are being held constant, 
the freestream Mach number of the flow drops with increasing gas constant. This decrease in Mach number causes 
lower post-shock temperatures which increases the induction time, shifting the behavior of the reaction front to an 
effective lower velocity compared to the baseline stoichiometric results. By the same logic a mixture that is fuel 
lean for a lighter-than-air fuel would have an effective higher velocity compared to the stoichiometric baseline. 
Note that the trends stated here would presumably be the opposite for a heavier than air fuel. 

V. Summary and Conclusions 

A quasi- ID formulation has been developed that provides equations describing axisymmetric flow valid in the 
stagnation region of a sphere. The equations used in this formulation were derived by applying a local similarity 
postulate to the axisymmetric Navier Stokes equations. In order to evaluate the validity of this formulation, a 
number of cases of spherical projectiles fired into hydrogen/air mixtures at hypersonic velocities were simulated 
inviscidly to compare against existing experimental and numerical studies. The first comparison showed a 
qualitative agreement between the quasi- ID simulation and selected experiments by Ruegg and Dorsey 4 that showed 
stable; regular, unsteady; and large-disturbance, unsteady behavior. In addition, large -disturbance results were 
compared to the numerical work of Matsuo and Fuj ii 2- which showed a qualitative agreement of the time evolution 
of the reaction and shock fronts. The second comparison took place over a range of projectile velocities 
encompassing the experiments conducted by Lehr' 1 . Numerical results by Sussman 23 were also included as an 
additional point of comparison. The velocity range was simulated densely by the quasi-ID code and at each point 
the frequency of the reaction front was recorded and compared to the experimental and higher-fidelity results. The 
general behavior of the reaction and shock fronts was captured over the velocity range in the present simulations and 
a reasonable agreement was obtained with recorded frequencies. The current simulations also showed trends not 
visible from the previous experimental and numerical results, warranting further investigation. Next, the simulated 
behaviors of the regular and large-disturbance instability regimes were compared against ID mechanism models 
present in the literature. It was found that all of the important flow features in the published models were captured 
by the quasi-ID formulation. Finally, the effect of equivalence ratio on flow behavior was evaluated. It was found 
that the reaction front was quite sensitive to inflow stoichiometry due to its effect on induction time. 

The present results are promising in establishing the quasi- ID formulation in rapidly evaluating premixed, 
shock-induced combustion flows, an important step toward understanding the behaviors of such flows over a large 
design space relevant to the operation of hypervelocity engines using PMSIC. Identifying regimes of predictably 
safe, steady behavior for conditions typically found at the entrance of the combustor is important in order to ensure 
feasible engine designs. In addition, understanding the mechanisms that cause large amplitude instabilities are 
helpful in avoiding or attempting to mitigate such behaviors. Overall, results from the quasi-ID analysis show 
promise as a valid starting point for analysts and designers in effectively conducting high-fidelity analyses and in the 
design of next-generation hypervelocity, airbreathing engines. 
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Figures 
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Figure 1. A schematic is shown for a hypervelocity vehicle using premixed, shock-induced combustion. 
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Figure 2. The problem setup for this study is shown. 
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Figure 3. a) Contours of water mass fraction on the stagnation streamline overlaid with contour lines of pressure as a 
function of non-dimensional time for Ruegg and Dorsey experiment Case 1. b) An annotated experimental shadowgraph 
captured by Ruegg and Dorsey for the same conditions is shown for comparison 4 . 
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Figure 4. a) Contours of water mass fraction on the stagnation streamline overlaid with contour lines of pressure as a 
function of non-dimensional time for Ruegg and Dorsey experiment Case 2. The drawn contact discontinuity is 
approximated based on separate density contours, b) An annotated experimental shadowgraph captured by Ruegg and 
Dorsey for the same conditions is shown for comparison 4 . 
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Figure 5. a) Contours of water mass fraction on the stagnation streamline overlaid with contour lines of pressure as a 
function of non-dimensional time for Ruegg and Dorsey experiment Case 3. b) An annotated experimental shadowgraph 
captured by Ruegg and Dorsey for the same conditions is shown for comparison 4 . 


16 

American Institute of Aeronautics and Astronautics 
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Figure 6. A detail view of density contours on the stagnation streamline as a function of time for Case 3 for the a) present 
quasi-ID simulation and b) high-fidelity, 2D simulation conducted by Matsuo and Fujii 22 (reproduced with permission). 
Note that the computed results are non-dimensional in time while the high-fidelity plot is in seconds. 
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Figure 7. Reaction front oscillation frequency as a function of projectile velocity for experimental and numerical 
investigations. 
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Figure 8. Contour lines of a) density on the stagnation streamline as a function of time for Ruegg and Dorsey experiment 
Case 2 and b) ID mechanism for the regular, unsteady regime by Matsuo 8 (reproduced with permission). 
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Figure 9. Contour lines of a) density on the stagnation streamline as a function of time for Ruegg and Dorsey experiment 
Case 2 and b) ID mechanism for the large-disturbance, unsteady regime by Matsuo and Fujii 22 (reproduced with 
permission). 
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Figure 10. A detail view of the a) stagnation streamline density contours from the present simulation (Case 3) compared 
to a b) ID sub-mechanism of the large-disturbance regime as proposed by Matsuo and Fujii 22 (reproduced with 
permission). 
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